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Abstract. The Bak-Tang-Wiesenfeld (BTW) sandpile model is a cellular automaton which has been inten- 
sively studied during the last years as a paradigm for self-organized criticality. In this paper, we reconsider 
a deterministic version of the BTW model introduced by Wiesenfeld, Theiler and McNamara, where sand 
grains are added always to one fixed site on the square lattice. Using the Abelian sandpile formalism we 
discuss the static properties of the system. We present numerical evidence that the deterministic model is 
only in the BTW universality class if the initial conditions and the geometric form of the boundaries do 
not respect the full symmetry of the square lattice. 



PACS. 64.60.Ht Dynamical critical phenomena - 05.65.+b Self-organized systems 
phenomena, random processes, noise, and Brownian motion 
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1 Introduction 

In equilibrium systems with short-ranged interactions and 
no broken continuous symmetry, correlation functions usu- 
ally decay exponentially with the distance. Exceptions 
from this behavior occur only in special cases, for example 
at critical points in a phase diagram. In this case the cor- 
responding correlation functions decay algebraically if the 
relevant system parameter is adjusted to special (critical) 
values (e.g. critical temperature, critical pressure, etc.). 

The concept of self-organized criticality (SOC) which 
was introduced by Bak, Tang and Wiesenfeld in 1987 MM 
attempts to explain the fact that in nature critical behav- 
ior is often observed, although nature cannot "fine-tune 
parameters" (see [||^ for an introduction and overview). 
The term "critical behavior" corresponds here to a power- 
law behavior of the probability distributions of certain 
physical quantities which characterize the system in both 
space and time |^,^. Typical examples for such quanti- 
ties are the size and life times of catastrophic events. The 
main idea is then that the critical state is an attractor of 
the dynamics. 

One of the paradigmatic systems which exhibit SOC 
is the Bak-Tang-Wiesenfeld sandpile model (BTW model) 
which was intensively investigated in the past. In the fol- 
lowing, we restrict our discussion to the BTW model on 
the two-dimensional square lattice. Concerning static prop- 
erties, analytical results exist which are mainly due to 
Dhar, who developed a formalism for Abelian sandpile 
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models [pi which allows to calculate exactly the height 
probabilities, height correlations, number of steady state 
configurations, etc. P,0;H;Pl. However, much less is known 
rigorously about the dynamical features, and estimates for 
the exponents of the probability distributions of avalanche 
quantities are only known from computer simulations (see 
for instance |l^,|ll|,|l2[|l^, lj,|l5[); further guesses for these 
exponents have been made via renormalization group ap- 
proaches (e.g. [[l6|p^ ). 

In this paper we consider a deterministic version of 
the BTW model which was introduced by Wiesenfeld et 
al. [|l8[ . In the original BTW model, sand grains are added 
at randomly chosen sites of the lattice. In the determinis- 
tic BTW model (DBTW) the seeding of sand is confined 
to one special site of the lattice. Computer simulations re- 
vealed [Esl that the DBTW model still displays criticality. 
Therefore, randomness in the location of the perturbations 
is not a necessary ingredient for SOC [Q. Furthermore, 
the authors concluded from their numerical analysis that 
the different versions of the BTW model could display 
different scaling behavior. 

However, this conclusion was obtained from an inves- 
tigation of the DBTW model for a small system size. As 
finite-size effects have been shown to affect the scaling 
behavior of the BTW model strongly [0,^3l ' ^^ reinves- 
tigate the case here and present a systematic finite-size 
analysis. We also present some new exact results for the 
DBTW model. 

The paper is organized as follows. In Sec. 0, we define 
the model and, using essentially the formalism developed 
by Dhar Q, study the static properties of the DBTW. 
In Sec. H, we present our results from computer simula- 
tions and then discuss the observed scaling behavior in 
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the context of the universahty hypothesis of Ben-Hur and 
Biham [|l9[. A summary closes the paper. 



2 Definition of the model and static 
properties 

The BTW model on a two-dimensional square lattice of 
size Lx L is defined as follows. To each site i {i = 1,2,..., 
L'^) an integer variable Zi (the height) is assigned. Starting 
with an empty lattice {zi = for all i) the addition of a 
grain of sand means to choose a site i at random and to 
increase Zi by 1. If Zi is equal to a fixed threshold value 
Zc, site i topples and distributes one grain of sand to each 
nearest neighbor, which can in turn trigger more toppling 
events. Thus, an avalanche of relaxation events may take 
place. For the sake of simplicity (and without loss of gen- 
erality) we set Zc = 4 throughout this work. The toppling 
rules can be formulated in terms oi a N x N toppling 
matrix A H, where N = L^. If site i topples, one has 



A,, 



(1) 



for all j = 1,2, ...,iV, where the toppling matrix satisfies 
the conditions An = 4, Z\y = —1 if i and j are nearest 
neighbors, and Aij = otherwise. Sand falling over the 
rim of the system is discarded. 

In the following we briefly describe the Abelian sand- 
pile formalism which was introduced by Dhar [pj and recall 
the major results. Dhar showed that the dynamics of the 
BTW model is well defined in the sense that the resulting 
stable configuration C = {zi} is always the same, regard- 
less of the order, in which critical sites are updated during 
an avalanche. One can define an operator a^ by its action 
onto a stable configuration C: UiC is the stable configu- 
ration which results from adding a particle at site i and 
relaxing the resulting configuration. Dhar showed that 



[ai,aj] = 0, 



(2) 



for all i,j. Therefore, the BTW model is called an Abelian 
sandpile model [|6|. 

One is interested in the stationary state of this cellu- 
lar automaton, i.e., one iterates the dynamical rules until 
all expectation values become time independent. Since the 
dynamics can be described as a Markovian process a sta- 
ble configuration can either be transient or recurrent. A 
recurrent configuration C can be defined by demanding 
that for every possible seeding site i a natural number 
mi{C) exists such that a"^ 'C = C holds. Thus, for a 
recurrent configuration C and for a natural number I one 
gets@ 

(3) 



aiC = oiaT^^'^^C = O^T'^'^MC, 



which shows that a[C is a recurrent configuration, too. 
Therefore, the set of recurrent configurations is closed un- 
der the action of the operators a; H . Since it is possible to 
define unique, inverse operators a^ , it follows that each 
recurrent configuration has the same probability to appear 



in the stationary state. One can further prove that the 
number of recurrent configurations is equal to det A |g] . 

A two point correlation function Gy can be defined 
in the following way: let Gy be the expectation value for 
the number of topplings in j , which are caused by adding 
a particle at site i. In the stationary state, the average 
number of particles which enter site j must be equal to 
the average of sand grains leaving site j: 



^ij'^jj ^ / ^ Gife(— Zlfej) -|- dij , 



(4) 



k^J 



which implies G = A ^ M. The analytic expression for G 
is well known: 
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with d = a7r/(L -|- 1), 6 = 67r/(L -I- 1) and where the sites 
i and j have the coordinates {xi,yi) and {xj,yj), respec- 
tively. Let s denote the size of an avalanche, i.e., the total 
number of topplings during that avalanche. Using Eq. (ra) 
it can be shown that the average number of topplings (s) 
is given by 



(=^)^ 



1 



L2 ^ 
i. 



L^{L- 



^ cot2 I cot^ I 



-l)'t:'sin2f + sin^| 

odd 



which scales for large L as 



{s)^L' 



(6) 



(7) 



Let us now turn to the deterministic BTW (DBTW) 
model as introduced by Wiesenfeld, Theiler, and McNa- 
mara |lq| . Here, sand is always added at a fixed input 
site io. Thus, the dynamics is fully deterministic and it is 
clear that in the configuration space (recurrent configura- 
tions) a DBTW model will settle down in an orbit with 
some period T, which in general will depend on jq. For ex- 
ample, the period of a DBTW model with a seeding site 
at the corner of the lattice is larger than the period of the 
center-seeded DBTW model, as it is intuitively clear and 
also known from computer simulations. An analytical ap- 
proach |2^] was used to evaluate T for the center-seeded 
DBTW model up to a system of A^ = 361 sites, where the 
authors found a period of length w 10^^ and extrapolated 
their results to reproduce the numerical estimation of [ |l8| 
T ^expjQ.llN). 

In ||l^ , several interesting features of the DBTW model 
could be derived by using the Abelian sandpile formalism: 

(i) The orbits have the same period (for a fixed input site 

ia), regardless of the initial conditions. 
(ii) The minimal stable configuration C* = {zi — Zc ~ 

1, for all i} is always on an orbit of the DBTW model. 
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Note that (ii) does not mean that a DBTW model will 
always reach C* at some point. For example, for the initial 
condition Zi — 1 for all i a system of 49 sites has an 
orbit without C* . Only for different initial conditions or 
a different system size (e.g. 25 sites with the same initial 
condition) C* is on the orbit. 

Let us now define Ni^j as the total number of topplings 
at site i within the period T of a DBTW model with input 
site i{). The time average of Ni^j is simply Ni^j/T. The 
total flow of particles during T into j has to equal the 
flow out of j and it follows Ni^j/T 



A;\ or 



^^J = G.„, = A^o,/T. 



(8) 



This means that every "row" i of the BTW correlation 
function Gij , which stands for the ensemble-averaged num- 
ber of topplings in j when seeding in i, represents the time 
average of topplings in j of a DBTW model with a fixed 
input site «o- The average avalanche size (s) for the BTW 
model can therefore be thought of as (sip) of the DBTW 
models averaged over all orbits and all possible seeding 
sites io — 1,2, ...,iV. 

Similar to the BTW model it is possible to calculate 
the average number of topplings (s) for the center-seeded 
DBTW model on a square lattice with length L [L odd). 
Using Eq. (||) we get with Xi^ = y^o = (L -I- 1)/2 
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with A = (2a + l)7r/(L -h 1) and B = {2a^ ^)'^I{L ^- 1), 
respectively. It is straightforward (though tedious) to eval- 
uate this expression for large L by using standard meth- 
ods. One obtains (s) ~ I? , analogous to the BTW model. 
It is also easy to get an expression for the time-averaged 
number of topplings Ni^i^/T of the input site zq: 



N,, 



T 



trj, 



(i + 1) 



^l? ^„«in2^ + 



a,fc=0 



sm" ^ + sin^ -j 



(10) 



which can be shown, again using elementary methods, to 
scale as — 7p — 



\nL for large L. 



3 Dynamics 

3.1 Characterization of the Avalanches 

For a graphical representation of the avalanches, it is con- 
venient to denote how many topplings n at each site have 
occurred during the avalanche |l^. Sites with the same 
number of topplings (during one avalanche) form "shells" . 
It is easy to check that for each site i inside such a shell. 



for which all four neighbors are also part of the shell, the 
height Zi before and after the avalanche remains the same. 
Figure |l| shows some examples which have been obtained 
in the stationary state of the central seeded DBTW. The 
shells seem to form compact sets with n monotonically 
decreasing from the central site towards the boundaries 
along the symmetry axis. Also, it seems that at the bound- 
ary of each n = const shell, n will change only by 1 (es- 
pecially at the boundary of the avalanche one finds n — 1 
most often). However, the latter statement is not always 
correct, there are rare exceptions. For instance, the upper 
left cluster in Fig. || has a boundary site which has toppled 
twice and the lower cluster has sites with n = 1 adjacent 
to a site with n = 3. We could only prove two properties 
of the avalanches. First, it has been shown for the BTW 
model, and the proof applies for the DBTW model also, 
that avalanches are always compact |1^,^. Second, one 
can show that n decreases monotonically from the cen- 
ter towards the boundaries along the symmetry axis by 
decomposing an avalanche into a series of waves of top- 
plings |^,23|. First one topples the center site and re- 
laxes all other sites which become unstable. This defines 
the first wave of topplings. After this, one allows the center 
site to relax again (if possible) which generates a second 
wave of topplings and so forth. It has been shown [^ 
that each avalanche produced by such a wave of topplings 
is compact. Furthermore, each site in a wave can topple 
only once. Thus, by superposition of the compact waves 
of topplings, n can only monotonically decrease from the 
center towards the boundaries along the symmetry axis. 



3.2 Scaling behavior of the avalanches 

In driven systems the boundary conditions can influence 
the stationary state |£J| . We show below that the scaling 
properties of the BTW model and its deterministic version 
are only the same if the latter has boundary conditions or 
an initial configuration which do not respect the square 
symmetry. 

We denote by s the total number of topplings which 
occurred during the lifetime t (in units of lattice sweeps) 
of an avalanche. The area a of such an avalanche is the 
number of distinct toppled sites. The outflow o is the to- 
tal number of sand grains which leave the system during 
an avalanche. The linear size of an avalanche is measured 
via the radius of gyration of an avalanche cluster. In the 
critical steady state the corresponding probability distri- 
butions should obey a power-law behavior 



Px{x) 



(11) 



characterized by the exponent t^ with x G {s, a, i, o, r}. As 
usual, we measure the avalanche distributions by counting 
the numbers of avalanches corresponding to a given area, 
duration, etc. and integrate these numbers over bins of 
increasing length (see for instance pSf ) . In our simulations 
successive bin lengths increase by a factor b — 1.2. In the 
case of the BTW model it is known that the probability 
distributions display logarithmic corrections 

v,-Tx const/ In L 



Px{x) 



(12) 
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the four possible values of s mod 4 (see Fig. O) . We have 
no explanation for this behavior. 

Let us briefly remark that a similar symmetry effect 
can be found in the outflow probability distribution Po{o). 
One gets different curves for all o which are divisible by 8 
and those which are not (see Fig. ^). To explain this, we 
note that the outflow o is divisible by 8 exactly when the 
sites in the middle of the boundary edges do not topple. 
This is extremely unlikely for large avalanches. On the 
other hand, for small avalanches, the constraint to topple 
sites in the middle of the boundary edges considerably 
reduces the number of possible avalanches. 

Since the exponents of the size and duration distri- 
bution, Ts and Tf, are not well defined here we renounce 
further investigations of the size and duration distribution 
in this section and focus our attention on the probability 
distributions of the avalanche area and radius which be- 
have as usual (see below). We measured the probability 
distributions Pa (a) and P^ (r) for various system sizes and 
obtained the corresponding exponents from a power-law 
fit of the straight portion of the curves. The values of both 
exponents are plotted in Fig. 0. In contrast to the BTW 
model [Eq.(|2|)] the avalanche exponents of the center- 
seeded DBTW model display no significant system size 
dependence. This allows us to apply the finite-size scaling 
analysis pq] 



P,ix,L) = L-^^ g^{xL-''^), 



(13) 



Fig. 1. Three different avalanches for the center-seeded 
DBTW model. Sites which toppled once are marked as A, twice 
toppled sites as B, etc. 



which are caused by finite-size effects [|Il|,|2|. A numerical 
determination of the avalanche exponents requires there- 
fore a careful analysis of finite-size effects. Using the func- 
tional form of the finite-size corrections it is possible to de- 
termine the exponents Tx directly, i.e. without any extrap- 
olation to the infinite system, and the best known values 
arer^ = 1.293±0.009, Ta == 1.33±0.011, r* = 1.480±0.011, 
andr^ = 1.665 ±0.013 H. 

We performed simulations of the center-seeded DBTW 
model for various system sizes L < 2049 and averaged 
all measurements over at least 2 x 10^ avalanches. Start- 
ing from an empty lattice, we added particles at the lat- 
tice center and applied the "burning algorithm" in order 
to check if the system has reached the steady state Q. 
Thus the initial conditions (as well as the boundary con- 
ditions) respect the square symmetry. Figure || shows the 
probability distribution Pt{t) of the avalanche duration. 
Surprisingly, it turns out that the above mentioned loga- 
rithmically averaging method is not suitable in our case, 
because avalanches of odd or even duration display a dif- 
ferent scaling behavior. The two branches of Pt{t) in Fig. 
{t°'^'^ and Tj°™") clearly have different slopes for the sys- 
tem size considered. The probability distribution of the 
avalanche size Ps{s) exhibits an even more complicated 
fine structure of four distinct branches corresponding to 



where the scaling exponents (3x and v^ are connected with 
the avalanche exponent t^ via the scaling equation f3x = 
TxVx [pql- This finite-size scaling ansatz works for the area 
and radius distribution and the corresponding data col- 
lapse for Pr{r) is plotted in Fig. p. The obtained values for 
the avalanche exponents agree with the results of the re- 
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Fig. 2. The probability distribution of the avalanche dura- 
tion Pt(t). Avalanches of an even or odd duration display a 
different scaling behavior. Thus, the usual logarithmic aver- 
aging of the distribution leads to useless results (solid line). 
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Fig. 3. The probability distribution of the avalanche dura- 
tion Ps{s). The distribution decomposes into four different 
branches corresponding to the four possible values of s mod 4. 



gression analysis (see Fig. ||) and we get Ta = 1.368±0.011 
andr^ = 1.752 ±0.027. 



As already mentioned in section 3.1, it is possible to 
show that the avalanches are compact. Thus, the area 
scales with the radius as 



(14) 



Then, the transformation law of probability distributions 
Pa{a)da — Pr{r)dr leads to the scaling relation 



This scaling relation is fulfilled within the error-bars, which 
confirms the accuracy of the determination of the avalanche 
exponents t^ and r^. Finally we mention that our obtained 
exponents are consistent with the values Ta = 11/8 and 
Tt = 7/4 and that these values obey the above scaling 
relation exactly. 



4 Discussion 

Since both avalanche exponents differ significantly from 
those of the BTW model we conclude that the center- 
seeded DBTW model does not belong to the BTW uni- 
versality class. It is worth to examine the different univer- 
sal behavior in detail since the universality hypothesis of 
Ben-Hur and Biham states that the only parameter which 
determines the scaling behavior (exponents) of a sandpile 
model is the so-called relaxation vector which describes 
how the sand grains of a critical site are distributed to 
Applying this concept of classifi- 



the next neighbors |19 



2 = 



1 



(15) 



cation one can identify three universality classes where the 
distribution is nondirected, nondirected on average, and 
directed. For instance, the BTW and the related Zhang 
model p^ belong to the universality class of nondirected 
models, whereas the Manna model [g^ is nondirected on 
average and therefore belongs to a different class. Sev- 
eral numerical investigations confirm these classification 
ansatz (see for instance [^,^|2^ and for recent investi- 
gations p^,^). 

According to this classification concept the center-seed- 
ed DBTW model and the BTW model should belong to 
the same universality class because both models are char- 
acterized by the same relaxation vector. In contrast to 
the BTW model the center-seeded DBTW model is deter- 
ministic and both the avalanches and the height config- 
urations display the square symmetry. Moving the input 
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Fig. 5. System size dependence of the avalanche exponents Ta 
Fig. 4. The probability distribution of the sand outflow Po{o) and r,.. The solid lines corresponds to the values obtained from 
for L = 129. The distribution decomposes into two different a finite-size scaling analysis and the dashed lines corresponds 
branches corresponding to the two possible values of o mod 8. to the values of the BTW model obtained in Q. 
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Fig. 6. The finite-size scaling analysis of the avalanche distri- 
bution Pr{r). Since the radius scales with the system size the 
exponent Vr should equal one. 



Fig. 8. The probability distribution Pais) of the two asym- 
metric DBTW models and the BTW model for various system 
sizes. For L < 513 the curves are shifted in the downward 
direction. 




Fig. 7. The probability distribution Pa{s) of the symmetric 
center-seeded DBTW model, two asymmetric DBTW models 
(in one case the square symmetry is broken by the boundaries 
and in the other case by asymmetric initial conditions), and 
the BTW model for L — 257. In the latter cases the curves are 
shifted in the downward direction. 



site Zo from the lattice center brakes the square symmetry 
but the dynamics of the system is still deterministic. The 
same effect is obtained if we use center-seeding but start 
with an asymmetric initial condition. We call these models 
the asymmetric DBTW models and our analysis revealed 
that they display the same scaling behavior as the usual 
BTW model. We plot the probability distribution Ps{s) 
for the considered models in Fig. u\ Except of deviations 
at the cut-off, the probability distributions of the asym- 
metric DBTW models agree with the corresponding curve 
of the BTW model and differ clearly from the distribution 



of the symmetric DBTW model. Figure ^ shows the prob- 
ability distribution for various system sizes. Again, apart 
from the cut-off behavior the curves of the BTW model 
and the asymmetric DBTW models are identical. This 
implies that the avalanche exponents of the asymmetric 
DBTW models display the same logarithmic corrections 
as the BTW model [Eq. (|l|)]. 

We conclude from our investigations that the asym- 
metric DBTW and the BTW model belong to the same 
universality class, whereas the center-seeded DBTW model 
with symmetric initial height configuration does not. Since 
the asymmetric DBTW model is still deterministic but 
lacks the square symmetry the different universal behav- 
ior of the center-seeded DBTW model is not caused by 
the deterministic dynamics but by the square symmetry 
of the system (in agreement with |lq| ). 

We conclude from our results that properties of the 
steady state such as symmetries or translational invari- 
ance can affect the universality class of sandpile models. 
This is confirmed by recently performed simulations [ pl| 
of a directed version of the Zhang model which exhibits 
a different scaling behavior than the exactly solved di- 
rected BTW model [B2|. According to the classification 
of Ben-Hur and Biham the directed Zhang model and the 
directed BTW model should belong to the same universal- 
ity class. But in contrast to the directed BTW model the 
height configuration of the directed Zhang model displays 
no translation invariance [ pl| . Similar to the center-seeded 
DBTW model one has to be careful to apply the univer- 
sality hypothesis of Ben-Hur and Biham. 

In summary we reconsidered a deterministic version of 
the Bak-Tang-Wiesenfeld sandpile model where the sand 
grains are added always to the central site of the lattice. 
Similar to the usual BTW model the Abelian sandpile for- 
malism allows to calculate some of the static properties 
of the system. Our numerical investigations show that the 
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deterministic central-seeded model with square symmetric 
initial conditions exhibits a different scaling behavior than 
the BTW model, in contrast to the deterministic model 
without square symmetry. 
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